clc; clear;

%% 1. 系统参数设置
params = struct(...
    'a1', 1.5, 'b1', 3, 'c1', 1, 'd1', 5, 'i1', 0.45, ...          % HR 神经元1
    'a2', 0.77, 'b2', 0.3, 'c2', 0.8, 'i2', 0, ...                   % FN 神经元
    'a3', 0.9, 'b3', 3, 'c3', 1, 'd3', 5, 'i3', 0.5, ...              % HR 神经元2
    'm12', 0.785, 'm21', 0.52, 'm23', 0.785, 'm32', 0.08);            % 耦合参数

% 固定其他初始条件
x2 = 3; x3 = 1; y2 = 1; y3 = 0; phi = 0.3;

eps_val = 4.87;

%% 2. 构造初始条件平面 (x1(0), y1(0))
x1_vals = linspace(-3, 3, 6);   % 调整网格数目以平衡精度与速度
y1_vals = linspace(-3, 3, 6);
[X1, Y1] = meshgrid(x1_vals, y1_vals);
[Nrows, Ncols] = size(X1);
N = numel(X1);

%% 3. 设置积分与 Lyapunov 计算参数
T_trans = 1000;    % 去除瞬态时间
T_le = 5000;       % LE 计算时间
dt_save = 5;      % 重正化时间间隔
RelTol = 1e-6;
AbsTol = 1e-9;
numLE = 1;        % 这里只需要计算最大 LE

% 预分配结果向量（按网格点展开）
LE_result = NaN(N,1);


parfor idx = 1:N
    % 根据索引获取对应初始条件
    xi = X1(idx);
    yi = Y1(idx);
    % 构造初始条件向量： [x1; y1; x2; y2; x3; y3; phi]
    X0 = [xi; yi; x2; y2; x3; y3; phi];
    LE_vals = LEs(eps_val, params, X0, T_trans, T_le, dt_save, RelTol, AbsTol, numLE);
    LE_result(idx) = LE_vals(1);
end

% 将结果重塑为二维矩阵
LE_map = reshape(LE_result, size(X1));

%% 5. 绘图
figure;
imagesc(x1_vals, y1_vals, LE_map > 0);  % 用逻辑值显示混沌 (LE > 0) 与周期 (LE <= 0)
axis xy; 
colormap([0 0 1; 1 1 0]);  % 蓝色表示周期，黄色表示混沌
colorbar;
xlabel('$x_1(0)$','Interpreter','latex','FontSize',12);
ylabel('$y_1(0)$','Interpreter','latex','FontSize',12);
title('Fig. 11(a): Attraction Basin (Blue: Periodic, Yellow: Chaotic)','FontSize',14);
set(gca, 'FontSize',12);
axis equal;
